This notebook produces the figures for the Springer Reference article "neuronal parameter space visualization". The code is GPL (see LICENSE.txt in top-level directory of the repository). If you use it to produce scientific results, please reference Schmuker, M (2013) Neuronal Parameter Space Visualization. Spinger Reference Computational Neuroscience, Springer, Berlin.
This notebook requires that you have PyNN and Brian installed. It was developed and tested with PyNN version 0.7.5 and Brian version 1.4.1, which you can get here: Brian: http://neuralensemble.org/trac/brian/downloader/download/release/19/1.4.1 PyNN: https://pypi.python.org/packages/source/P/PyNN/PyNN-0.7.5.tar.gz
The most recent versions can be found on http://neuralensemble.org.
In [11]:
import sys
import os
sys.path.append(os.getcwd())
from paramspacevisu import network
First we define the network class. Two neuronal populations, modeled as poisson spike sources, projecting onto a third population, the target population.
In [ ]:
default_params = { #default parameters for the network
"wi": 0.01, # inhibitory weight
"we": 0.005, # excitatory weight
"ri": 30., # inhibitory rate
"re": 30., # excitatory rate
"p_conn_et": 0.1, # connection probability exc-> target
"p_conn_it": 0.1, # connection probability inh-> target
"num_e":200, # number of excitatory neurons
"num_i":50, # number of inhibitory neurons
"num_t":30} # number of target neurons. Reduce this for faster simulation with very similar results.
class ThreePops(object):
"""
Three populations of neurons: Excitatory, inhibitory, and target.
200 exc (poisson), 50 inh (poisson), 50 target (IF).
convergence: random 10 % (exc to target and inh to target)
weights: w_inh = 4*w_exc
"""
def make_network(self, param_dict=default_params):
"""
Construct the network according to the parameters specified.
"""
pynn.setup()
self.exc = pynn.Population(param_dict['num_e'],
pynn.SpikeSourcePoisson,
cellparams={'rate':param_dict['re']})
self.inh = pynn.Population(param_dict['num_i'],
pynn.SpikeSourcePoisson,
cellparams={'rate':param_dict['ri']})
self.target = pynn.Population(param_dict['num_t'],
pynn.IF_cond_exp)
self.target.record(to_file=False)
connector_et = pynn.FixedProbabilityConnector(
p_connect=param_dict["p_conn_et"],
weights=param_dict["we"])
connector_ei = pynn.FixedProbabilityConnector(
p_connect=param_dict["p_conn_it"],
weights=param_dict["wi"])
self.prj_et = pynn.Projection(self.exc, self.target,
method=connector_et)
self.prj_it = pynn.Projection(self.inh, self.target,
method=connector_ei, target='inhibitory')
def update_weights_and_rates(self, param_dict):
"""
Update the weights and firing rates in the network according to the
parameter dict.
"""
self.exc.set('rate', param_dict['re'])
self.inh.set('rate', param_dict['ri'])
self.prj_et.setWeights(param_dict['we'])
self.prj_it.setWeights(param_dict['wi'])
def run_network(self, duration=1000.):
"""
Run the network for the specified duration. Returns spikes of target
population in GDF format.
"""
pynn.run(duration)
spikes = self.target.getSpikes()
return spikes
If you just want to play around with the visualization of parameters, the next four cells produce some examples for 2D slices through the parameter space. The code to reproduce the figures in the reference article can be found after these examples.
In [5]:
#make one panel: rate variation
we = 0.005
wi = 0.01
num_spikes_1p = numpy.zeros((5,5))
param_dicts_1p = numpy.zeros((5,5), dtype=object)
net = network.ThreePops()
net.make_network(param_dict=default_params)
lastnumspikes = 0
for ire,re in enumerate(numpy.linspace(20.,40.,5)):
for iri,ri in enumerate(numpy.linspace(20.,40.,5)):
print "re:%.1f, ri:%.1f, we:%.4f, wi:%.4f"%(re,ri,we,wi)
params = default_params.copy()
params.update({'re':re, 'ri':ri, 'we':we, 'wi':wi})
param_dicts_1p[ire,iri] = params
net.update_weights_and_rates(params)
spikes = net.run_network()
currentnumspikes = len(spikes)
num_spikes_1p[ire,iri] = (currentnumspikes - lastnumspikes)/params['num_t']
lastnumspikes = currentnumspikes
In [6]:
imshow(num_spikes_1p, interpolation='nearest', origin='lower')
colorbar()
Out[6]:
In [7]:
#make one panel: weight variation
re = 30.
ri = 30.
num_spikes_1pw = numpy.zeros((5,5))
param_dicts_1pw = numpy.zeros((5,5), dtype=object)
net = network.ThreePops()
net.make_network(param_dict=default_params)
lastnumspikes = 0
for iwe,we in enumerate(numpy.linspace(0.003,0.007,5)):
for iwi,wi in enumerate(numpy.linspace(0.008,0.012,5)):
params = default_params.copy()
params.update({'re':re, 'ri':ri, 'we':we, 'wi':wi})
param_dicts_1pw[iwe,iwi] = params
net.update_weights_and_rates(params)
spikes = net.run_network()
currentnumspikes = len(spikes)
num_spikes_1pw[iwe,iwi] = (currentnumspikes - lastnumspikes)/params['num_t']
lastnumspikes = currentnumspikes
In [8]:
imshow(num_spikes_1pw, interpolation='nearest', origin='lower')
colorbar()
Out[8]:
This cell samples the parameter space of the network systematically along the four parameters $r_{exc}$, $r_{inh}$, $w_{exc}$, and $w_{inh}$. Each parameter is sampled at five equal-spaced intervals.
Note: The simulations take a long time to run. The progress is displayed by the parameter combinations that have already been evaluated. You can skip the simulation and load the resulting data directly in the next cell below.
In [27]:
# first panel: vary re and ri between 10 and 50 in five steps
num_tn = 50.
num_spikes = numpy.zeros((25,25))
param_dicts = numpy.zeros((25,25), dtype=object)
net = network.ThreePops()
net.make_network(param_dict=default_params)
lastnumspikes = 0
for ire,re in enumerate(numpy.linspace(20.,40.,5)):
for iri,ri in enumerate(numpy.linspace(20.,40.,5)):
for iwe,we in enumerate(numpy.linspace(0.003,0.007,5)):
for iwi,wi in enumerate(numpy.linspace(0.008,0.012,5)):
print "re:%.0f, ri:%.0f"%(re,ri)
params = default_params.copy()
params.update({'re':re, 'ri':ri, 'we':we, 'wi':wi})
index = (iwe*5 + ire, iwi*5+iri)
param_dicts[index] = params
net.update_weights_and_rates(params)
spikes = net.run_network()
currentnumspikes = len(spikes)
num_spikes[index] = (currentnumspikes - lastnumspikes)/num_tn
lastnumspikes = currentnumspikes
import cPickle
cPickle.dump(num_spikes, open('numspikes.cPickle','w'))
cPickle.dump(param_dicts, open('paramdicts.cPickle','w'))
Simply load the data if you don't want to wait for the simulations by executing this cell.
In [2]:
import os
import cPickle
datapath = os.path.realpath(os.path.join(os.getcwd(), os.path.pardir))
num_spikes = cPickle.load(open(os.path.join(datapath, 'numspikes.cPickle'),'r'))
param_dicts = cPickle.load(open(os.path.join(datapath, 'paramdicts.cPickle'),'r'))
We produce the first 2D-slice through the parameter space by picking all entries in the simulation where the firing rates $r_{exc}$ and $r_{inh}$ change while the weights $w_{exc}$ and $w_{inh}$ remain fixed.
In [42]:
vmax = numpy.max(num_spikes) # needed for consistent colorbar scaling
slice1_ind = [(x,y) for x in range(25) for y in range(25)
if (abs(param_dicts[x,y]['we']-0.005) < 0.0001)
and (abs(param_dicts[x,y]['wi']-0.010) < 0.00001)]
#print slice1_ind
#len(slice1_ind)
rates_slice1 = numpy.array([num_spikes[sl1] for sl1 in slice1_ind]).reshape((5,5))
#print rates_slice1
imshow(rates_slice1, interpolation='nearest', origin='lower', vmax=vmax)
ax = pylab.gca()
ax.set_xticklabels(['','20','25','30','35','40'])
ax.set_xlabel('$r_{inh}$ [spikes/s]')
ax.set_yticklabels(['','20','25','30','35','40'])
ax.set_ylabel('$r_{exc}$ [spikes/s]')
#ax.set_title('$w_{exc}=0.005$, $w_{inh}=0.010$')
cb = colorbar()
cb.set_label('spikes/s')
f = pylab.gcf()
f.set_size_inches(3,2.5)
f.tight_layout()
f.savefig(os.path.join(os.path.pardir, 'images', 'panelB.png'), dpi=600)
Similarly we proceed for a 2D-slice that shows the behavior of the network for constant firing rates with changing weights
In [43]:
slice2_ind = [(x,y) for x in range(25) for y in range(25) if (param_dicts[x,y]['re']==30.) and (param_dicts[x,y]['ri']==30.)]
#len(slice2_ind)
rates_slice2 = numpy.array([num_spikes[sl2] for sl2 in slice2_ind]).reshape((5,5))
#print slice2_ind
#print numpy.array(slice2_ind,dtype=object).reshape((5,5,2))
#print rates_slice2
imshow(rates_slice2, interpolation='nearest', origin='lower', vmax=vmax)
ax = pylab.gca()
ax.set_xticklabels(['','8','9','10','11','12'])
ax.set_xlabel('$w_{inh}$ [fS]')
ax.set_yticklabels(['','3','4','5','6','7'])
ax.set_ylabel('$w_{exc}$ [fS]')
#ax.set_title('$r_{exc}=30$, $r_{inh}=30$')
#cb = colorbar()
#cb.set_label('spikes/s')
f = pylab.gcf()
f.set_size_inches(2.45,2.45)
f.tight_layout()
f.savefig(os.path.join(os.path.pardir, 'images', 'panelC.png'), dpi=600)
Now we show the entire dimensional stack for all four parameters at once (unfortunately without axis annotations).
In [44]:
imshow(num_spikes, interpolation='nearest', origin='lower', vmax=vmax)
#cb = colorbar()
#cb.set_label('spikes/s')
#x-axis: fast: ri 20->40, slow wi 0.008->0.012
#y-axis: fast: re 20->40, slow we 0.003->0.007
ax = pylab.gca()
ax.set_xticks([])
ax.set_yticks([])
f = pylab.gcf()
f.set_size_inches(5,5)
f.tight_layout()
f.savefig(os.path.join(os.path.pardir, 'images', 'panelD.png'), dpi=600)